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The random walk problem is studied in two and three di- 
mensions in the presence of a random distribution of static 
traps. An efficient Monte Carlo method, based on a mapping 
onto a polymer model, is used to measure the survival proba- 
bility P(c, t) as a function of the trap concentration c and the 
time t. Theoretical arguments are presented, based on earlier 
work of Donsker and Varadhan and of Rosenstock, why in two 
dimensions one expects a data collapse if — ln[P(c, t)]/ \n(t) 
is plotted as a function of \/A7/ln(t) (with A = — ln(l — c)), 
whereas in three dimensions one expects a data collapse if 
-t~ 1/3 ln[P(c, t)] is plotted as a function of t 2/3 X. These ar- 
guments are supported by the Monte Carlo results. Both data 
collapses show a clear crossover from the early-time Rosen- 
stock behavior to Donsker- Varadhan behavior at long times. 

The problem of trapping in a random medium is funda- 
mental for understanding many important physical pro- 
cesses in disordered systems. Practical applications in- 
clude some dynamical processes in disordered media, ki- 
netics of reactions, electron-hole recombinations in ran- 
dom and amorphous solids, exciton trapping and anni- 
hilation, etc. In the lattice version of this problem, 
traps are randomly distributed on a lattice of dimension 
d. Usually, all correlations between traps are ignored, 
largely motivated by the fact that in most relevant phys- 
ical situations the concentration of traps is small. Ini- 
tially at time t = 0, a large population of random walk- 
ers is uniformly distributed over the lattice. Each walker 
hops from one lattice site to a randomly chosen nearest- 
neighbor site, at a rate of one hop per time unit. When 
a walker meets a trap, it dies. The time dependence of 
the survival probability of the walkers in the presence of 
traps is an interesting problem, which has been analyzed 
both mathematically and numerically in the literature 
|p]-|l5|. The problem is related to the properties of the 
density of states of the electrons on a lattice with ran- 
domly distributed impurities, which has been discussed 
in Refs. 

Theoretical treatments of the trapping problem usually 
start by considering the configuration averaged survival 
probability P(c, t), where c is the concentration of static 
traps. If we assign a label i = 1 . . . N to all different ran- 
dom walks of t steps and introduce n(i,t) as the number 
of different sites visited by the walk labeled i we can write 
the survival probability as 



P(c,i)H(l-c)" (M) ) 



(1) 



Rosenstock (RS) M expression 

(1 _ C )WM)>. This leads 
to exponential dependence on t for d > 3 and exponential 
dependence on \fi in d = 1. The short-time limit in two 
dimensions will be discussed further on in this paper. 

The behavior of the survival probability at long times 
has been treated rigorously by Donsker and Varadhan [Q ; 
the limit t — > oo is known in the literature as Donsker- 
Varadhan (DV) limit. In the DV limit, the survival prob- 
ability P(c, t) in d dimensions does not decay exponen- 
tially, as simple intuition would suggest, but in a rather 
more complicated way: 



ln[P(c,t)] 



-A\ d + 2 t d + 2 . 



(2) 



where A is a constant depending on dimension as well 
as on the characteristics of the lattice and the random 
walk process, and A = — ln(l — c), with c the concentra- 
tion of static traps. An interpretation of this apparently 
unusual behavior was given by Grassberger and Procac- 
cia ||: they assign it to the existence of very rare large 
trap-free regions where walkers can survive for a long 
time. With increasing time ever larger trap free regions 
become dominant; the probability of finding such regions 
decreases exponentially with their d-dimensional volume, 
but the decay rate of particles moving within such a re- 
gion is inversely proportional to the square of its diame- 
ter. The optimal choice of this diameter gives rise to the 
stretched exponential behavior of Eq. (||) . 

One important problem with the DV result is that it 
does not say anything about how large t should be, in 
order to observe this behavior, nor what should be the 
corresponding value of the survival probability. The is- 
sue of the onset of the DV behavior has been studied 
numerically by a variety of algorithms, resulting in esti- 
mates which vary widely, apparently depending on the 
kind of algorithms used in the calculation. Earlier re- 
sults of Klafter et al. JtT], based on their numerical sim- 
ulations, suggest that P(c,t) should be less than 10~ 21 
for d — 2 to observe this behavior, whereas Fixman || 
obtained a value of P(c,t) < 10~ 67 in d = 3. Among 



since the probability that a site is not a trap is equal to 
(1— c). For short times P(c, t) is well approximated by the 



others, Havlin et al. ||9|] used an exact enumeration ap- 
proach and obtained P(c,t) < 10~ 13 for both two and 
three dimensions. The scanning method of Meirovitch 
ifiof supported the claims of Havlin et al. None of the 
above studies have shown conclusively the crossover from 
the early-time RS behavior to the asymptotic DV be- 
havior, except perhaps at very high trap densities. The 
most recent progress comes from Gallos et al. Jl4|], who 
have studied the problem numerically in d = 2, using 
the idea of self-interacting random walks coupled with a 
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"slithering snake" algorithm. They plotted the quantity 
— 1~ 01 ln[P(c, t)] as a function of t 8 X and obtained a 
fairly satisfactory collapse, within their numerical preci- 
sion; the theoretical basis for this kind of data collapse 
is however unclear. 

Our algorithm for the symmetric random walk on a 
simple cubic or quadratic lattice, the case studied mostly 
before, is also based on the identification of each random 
walk of t steps with a polymer configuration of length t. 
A direct Monte Carlo procedure to estimate the survival 
probability would then be to generate a large number 
of random walks, and average the quantity (1 — c) n ^ l,t \ 
This procedure would however give poor statistics with 
increasing t, since the result would be dominated by the 
one or two random walks that are the most compact. 
More accurate Monte Carlo results can be obtained by 
biasing the random walks generated towards more com- 
pact ones. Here, the identification of random walks with 
polymers comes into play. To each polymer configuration 
i we attribute an energy, given by the Hamiltonian 



H(i,t) = Xn(i,t). 



(3) 



If we then sample polymer configurations at infinite tem- 
perature, the average energy is equal to \n(i,t), which in 
turn equals ln[P(c, t)]. Numerically more accurate results 
are obtained, however, by generating polymer configura- 
tions at some inverse temperature (3. The likelihood for a 
configuration i to be generated is then proportional to its 
Boltzmann factor exp(—(3H(i,t)) and as a consequence 
there is a bias towards compact configurations. Also from 
this ensemble, one can compute the survival probability 
using 



P(c,t) = 



(cxp{+^H(i,t)){l-c) n ^) 
{exp(+pH(i,t)))p 



(4) 



where the averaging is performed over the ensemble gen- 
erated at inverse temperature (3. Note that at infinite 
temperature, i.e. (3 = 0, there is no bias and we retrieve 
Eq.(|lj). Interestingly, at inverse temperature (3 — 1 the 
numerator equals unity and the inverse of the survival 
probability equals the average of the inverse Boltzmann 
weights. It turns out that the most accurate results are 
obtained at some inverse temperature between and 1; 
all our results are obtained at (3 — 0.75. 

We generate a Markov chain of polymer configurations 
at inverse temperature (3 by proposing randomly cho- 
sen single-monomer moves, and accepting each proposed 
move with the Metropolis acceptance probability 



min [1, exp(-PAE)] , 



(5) 



where AE is the change in energy, defined by the Hamil- 
tonian Eq. (J^j . The types of single-monomer moves that 
we propose are depicted in Fig. [I]. In a naive implemen- 
tation, the computational effort per move would scale 
linearly with polymer length because of the necessity to 
walk along the polymer to compute the change in energy. 



However, if one stores in each lattice site the number of 
times the polymer visits that site, then the energy dif- 
ference can be computed easily: the energy increases by 
A if the monomer moves to a site which is not visited 
before, decreases by A if the monomer leaves a site that 
it visits only once, and stays unaltered if either both or 
neither of these two situations hold. This administration 
is also easily updated: the new site is visited one more 
time, while the old site is visited one time less. 
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FIG. 1. Monte Carlo moves used in the simulation: ex- 
change of a pair of consecutive steps (left and middle), or 
replacement of the first or last step (right). 

We now discuss numerical results obtained using our 
algorithm outlined above. First, we concentrate on the 
case of two dimensions. In Fig. [| we have plotted the log- 
arithm of the survival probability P(c, t) for concentra- 
tions c = 0.01, 0.02, . . . , 0.09; 0.1, 0.2, . . . , 0.9, and times 
(polymer lengths) up to t = 825. In these simulations, 
we start with a random walk, which is then simulated 
over 10 9 Monte Carlo steps per monomer, at an inverse 
temperature of (3 = 0.75. Data obtained during the first 
50000t attempted Monte Carlo steps per monomer have 
been discarded (thermalization) . 
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FIG. 2. Monte Carlo results for the survival probability 
as a function of time, for several trap concentrations on a 
two-dimensional lattice: — ln[P(c, t)] is plotted as a function 
of time t. 
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We now obtain a theoretically justified data collapse 
based on the identification of a common scaling variable 
for the DV and RS regimes. At short times and small 
concentrations, we expect that the polymer configura- 
tion closely resembles a random walk. In the case of two 
dimensions, this means that the number of different sites 
visited scales as t/ln(t) jL8|. Consequently, we deduce 
that the RS behavior in two dimensions is properly de- 
scribed by 



- ln[P(c,t)] ~ At/lnO). 



(6) 



At long times and large concentrations we expect to ob- 
serve the DV behavior: 



ln[P(c,f)]~(At) 



1/2 



(7) 



To obtain the proper scaling variable, we equate the right 
terms in equations (||) and (0), and obtain VXt/ ln(t) = 1; 
we can use this as a scaling variable, and rewrite equa- 
tions (H) and (0) as: 



-ln[P(c,t)] 



In(t) 

MPM)] 
ln(t) 



ln(t) 
Xt 



Hty 



(8) 



Thus we expect that if — ln[P(c, i)]/ ln(t) is plotted as 
a function of vAi/hi(t), the data for all values of trap 
concentration c and time t collapse onto a single curve, 
with an effective exponent that crosses over from 2 to 1. 

In Fig. |3| we have performed this plot, using the same 
data as in Fig. 0. Clearly, the data collapse is convincing 
over the whole range of parameters used in our simula- 
tion. Also, one can clearly identify the RS regime with 
a slope of 2 and the DV regime with a slope of I, both 
indicated with solid lines. The numerical estimate of the 
point where these lines cross is — ln[P(c, t)]/ ln(i) = 3.5 
and vAi/ln(t) = 1.13. The survival probability where 
the DV regime starts is then given by — ln[P(c, t)] = 
S.lVXt, for a suitable choice of c and t. The fact that 
the survival probability at the start of the DV regime is 
not a constant explains the wide range of reported values 
for this quantity. The apparent dependence of the results 
on the simulation methods can be understood because 
certain simulation methods are especially suited for high 
c-values whereas other methods require fairly low values 
of c. From our results one may conclude that the DV 
regime can already start at a significant survival proba- 
bility, depending on c; it should therefore be observable 
experimentally. 

To illustrate this in detail we plotted in Fig. || the 
cross-over values of both time and survival probability as 
functions of the trap concentration. 
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FIG. 3. Collapse of the two-dimensional data: 
— ln[P(c, t)]/ ln(t) is plotted as a function of yf~Xi/\n{i) i n a 
double-logarithmic plot. The solid lines are fits to the data, 
with slopes 2 and 1. They cross at the point (1.13, 3.5). 




FIG. 4. Location of the crossover probability P(c, i) as a 
function of trap concentration c, in two dimensions. The inset 
shows the crossover time t as a function of c. 

In dimensions different from d — 2 it is equally pos- 
sible to define a common scaling variable for the RS 
and DV regimes, that ought to give rise to a data col- 
lapse. For d < 2 one should consider ln[P(c, t)] as a func- 
tion of (\t) d / 2 and for d > 2 collapse should occur for 
_t(2/<*-l)ln[p( C) t)] as a function of t 2 ' d \. From these 
expressions one first of all sees very clearly that d = 2 
acts as a cross-over value for the dimensionality, as sug- 
gested already by the logarithmic terms in the 2d scaling 
variables. 

In particular, in three dimensions, where the number 
of different sites visited increases linearly with time t, the 
RS and DV regimes can be written as: 



,-1/3 



hx[P(c,t)} =< 2/3 A; 
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-1/3 



ln[P(c,t)} = 



t 2 / 3 x 



2/5 



(9) 



To show the resulting data collapse in three dimen- 
sions, we plotted in Fig. | -t" x /3 l n [P(c, t)} func- 
tion of t 2 / 3 \. All data have been obtained again with 
the polymer algorithm described above. As in two di- 
mensions, they can be explained by a crossover from the 
RS regime, where the curve has a slope of 1, to the DV 
regime, where the curve has a slope of 2/5. 



In summary, we have outlined an efficient algorithm 
to investigate diffusion with random traps in an arbi- 
trary concentration range. The simulation results for all 
trap concentrations and times could be explained by a 
crossover from the Rosenstock behavior at short times to 
the Donsker-Varadhan behavior at long times. We have 
identified the location of the crossover point, which in 
contrast to earlier belief, does not simply yield a specific 
value for the survival probability. 
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FIG. 5. Collapse of the three-dimensional data: 
— i -1 / 3 ln[P(c, t)] is plotted as a function of t 2//3 A in a dou- 
ble-logarithmic plot. The solid lines are fits to the data, with 
slopes 1 and 2/5. They cross at the point (8.5, 6.3). 

In Fig. we also plotted the cross-over time and sur- 
vival probability as functions of c for three dimensions. 
As one might expect, the concentration dependence is 
even much stronger than in two dimensions. 
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FIG. 6. Location of the crossover probability P(c, i) as a 
function of trap concentration c, in three dimensions. The 
inset shows the crossover time t as a function of c. 



[1] S. Havlin, D. Ben-Avraham, Adv. Phys. 36, 696 (1987) 

and references therein. 
[2] J. R. Beeler, Phys. Rev. A 134, 1396 (1964). 
[3] J. Haus and K. W. Kehr, Phys. Rep. 150, 263 (1987). 
[4] M. Donsker and S. R. S. Varadhan, Commun. Pure 

Appl. Math. 28, 525 (1975); ibid 32, 721 (1979). 
[5] P. Grassberger and I. Procaccia, J. Chem. Phys. 77, 6281 

(1982). 

[6] H. B. Rosenstock, J. Math. Phys. 11, 487 (1970). 

[7] J. Klafter, G. Zumofen and A. Blumen, J. Phys. (Paris) 

45, L49 (1984). 
[8] M. Fixman, Phys. Rev. Lett. 52, 791 (1984). 
[9] S. Havlin et. al, Phys. Rev. Lett. 53, 407 (1984). 
[10] H. Meirovitch, Phys. Rev. A 32, 3699 (1985). 
[11] Th. M. Nieuwenhuizen. Phys. Rev. Lett. 62, 357 (1989). 
[12] T. C. Lubensky, Phys. Rev. A 30, 2657 (1984); A. B. Har- 
ris and T. C. Lubensky, Phys. Rev. B 23, 2640 (1981). 
[13] J. K. Anlauf, Phys. Rev. Lett. 52, 1845 (1984). 
[14] L. K. Gallos, Panos Argyrakis and K. W. Kehr, 

Phys. Rev. E 63, 21104 (2001). 
[15] F. den Hollander and G. H. Weiss, in Contemporary Prob- 
lems in Statistical Physics, edited by G. H. Weiss (SIAM, 
Philadelphia, 1994). 
[16] I. M. Lifshitz, Usp. Fiz. Nauk. 83, 617 (1964); 

Sov. Phys. Usp. 7, 549 (1965). 
[17] R. Friedburg and J. M. Luttinger, Phys. Rev. B 12, 4460 
(1976). 

[18] E. W. Montroll and G. H. Weiss, J. Math. Phys. 6, 167 
(1965); J. W. Haus and K. W. Kehr, Phys. Rep. 150, 
263 (1987). 

[19] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, 
A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 
(1953). 



4 



